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Abstract. - Conserved growth models that exhibit a nonlinear instability in which the height 
(depth) of isolated pillars (grooves) grows in time are studied by numerical integration and 
stochastic simulation. When this instability is controlled by the introduction of an infinite 
series of higher-order nonlinear terms, these models exhibit, as function of a control parameter, 
a non-equilibrium phase transition between a kinetically rough phase with self-affine scaling 
and a phase that exhibits mound formation, slope selection and power-law coarsening. 



The nonequilibrium kinetics of the growth of films by the deposition of atoms on a substrate 
is of considerable experimental and theoretical interest 0111 . While the process of kinetic 
roughening pP leading to a self-afRnc interface profile has been extensively studied, there has 
been much recent interest ,2^S'A. ^ in a different mode of surface growth, involving the formation 
and coarsening of "mounds" (pyramid- like structures). The system is said to exhibit slope 
selection if the typical slope of the sides of the mounds remains constant during the coarsening 
process. Traditionally, the formation of mounds has been attributed to the presence of an 
Ehrlich-Schwoebel (ES) step-edge barrier that hinders the downward motion of atoms 
across the edge of a step. The destabilizing effect of the resulting "uphill" surface current 
is usually modeled in continuum growth equations [3|S| as a linear instability arising from a 
Laplacian of the height variable with a negative coefficient. 

In this Letter, we show that mound formation and power-law coarsening with slope selec- 
tion occurs in a class of well-known, conserved surface growth models as a result of a nonlinear 
instability which leads to a dynamical phase transition between kinetically rough and mounded 
morphologies. We consider the conserved, fourth-order, nonlinear growth equation proposed 
by Lai and Das Sarma (Oj and by Villain 

dh'{r, t')/dt' = -vV^h' + A'V^I V/if + ry'(r, t'), (1) 
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where h'{r,t') represents the height variable at the point r at time t' , V and represent, 
respectively, the spatial derivative and Laplacian operators in d dimensions (the dimension 
of the substrate), and 77' is a Gaussian, delta-correlated random noise. This equation is 
believed (2ii9nl0) to provide a correct description of the scaling behavior of kinetically rough 
surfaces of films grown by molecular beam epitaxy. Our results are based on numerical 
integration of this equation in one dimension, using a simple Euler scheme Choosing 
appropriate units of height, length and time, and discretizing both space and time, Eq.QJ is 
written as 

hi{t + At) - hi{t) = AtV^[-V^hi{t) + X\Vhi{t)f] 

+ VMf],{t), (2) 

where hi (t) represents the dimensionless height variable at the lattice point i at dimensionless 
time t, V and are lattice versions m] of the derivative and Laplacian operators, and r]i{t) 
is a random variable with zero average and variance equal to unity. These equations, with an 
appropriate choice of At, are used to numerically follow the time evolution of the interface. 
We have also studied an atomistic model ^3 in which the height variables {hi} are integers. 
This model is defined by the following deposition rule. First, a site (say i) is chosen at random. 
Then the quantity 

K,{{h,})^-V^h, + \\Vh,\^ (3) 

is calculated for the site i and all its nearest neighbors. Then, a particle is added to the site 
that has the smallest value of K among the site i and its nearest neighbors. In the case of 
a tie for the smallest value, the site i is chosen if it is involved in the tie; otherwise, one of 
the sites involved in the tie is chosen randomly. The number of deposited layers provides a 
measure of time in this model. 

It was found earlier that both these models exhibit a nonlinear instability in which 
isolated structures (pillars for A > 0, grooves for A < 0) grow rapidly if their height (depth) 
exceeds a critical value. This instability can be controlled T°P by replacing |V/iip in Eqns.lO 
and l|2Jl by /(|V/iip) where the nonlinear function f{x) = (1 — e~'^^)/c, c > being a 
control parameter. We call the resulting models "model I" and "model 11" , respectively. This 
replacement, which amounts to the introduction of an infinite series of higher-order nonlinear 
terms, is physically meaningful. Politi and Villain ^ have shown that the nonequilibrium 
surface current that leads to the V^|V/i'p term in Eq.QJ should be proportional to V|V/i'p 
when |V/i'| is small, and should go to zero when |V/i'| is large. The introduction of the 
"control function" /(|V/iip) satisfies this physical requirement. 

The time evolution of the height variables in model I is, thus, given by 

h.,{t + At) - h,{t) ^ AtV^[~V^h,{t) 

+ A(l-e-^l^'''(*)l')/c]+VAt?7,(i). (4) 

It was shown in Ref. that if c is sufficiently large, then the instability is absent and both 
models exhibit the scaling behavior expected in kinetic roughening. In the present work, 
we show that as the value of c is decreased, these models exhibit a first-order dynamical 
phase transition |13| to a mounded morphology at a critical value of c. The mounded phase 
exhibits power-law coarsening (interface width W ~ ): while the slope of the mounds 
remains constant. We present results for the phase diagram of these models in the (A, c) plane 
and describe the results of a stability analysis that provides an understanding of the observed 
behavior. We also show that this mechanism of mound formation is qualitatively different from 
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the conventional ES mechanism. Our results are quite general in that the behavior found in 
our models does not depend crucially on the form of the function f{x): any monotonic function 
that is linear for small x and saturates for large x leads to similar results. In particular, we 
have found very similar behavior using the form, f{x) = x/ [1 + ex), suggested by Politi and 
Villain 0. 

Our results are obtained for systems of different sizes (40 < L < 1000 - we do not find 
any significant dependence of the results on L) with periodic boundary conditions. In most 
of our studies of model I, we used At = 0.01. We have checked that very similar results are 
obtained for smaller values of At. If the control parameter c is sufficiently large, then the 
nonlinear instability is completely suppressed and the models exhibits the usual dynamical 
scaling behavior with the expected exponent values, (3 ~ 1/3, the dynamical exponent 
z ~ 3, and the exponent a = /3z ~ 1. As the value of c is decreased with A held constant, 
the instability makes its appearance: the height ho of an isolated pillar (for A > 0) increases 
in time if hmm{\,c) < ho < hrnax{^,c). The value of hmin is nearly independent of c, while 
hmax increases as c is decreased. If c is sufficiently large, hmax is small and the instability 
does not affect the scaling behavior of global quantities such as W. As c is decreased further, 
hmax becomes large, and when isolated pillars with ho > h„iin are created at an initially flat 
interface through random fluctuations, the rapid growth of such pillars to height h^ax leads 
to a sharp upward departure from the power-law scaling of W with time t. The time at which 
this departure occurs varies from run to run. Typical results obtained for model I with A — 4.0 
and c = 0.02 are shown in the inset of Fig^ 

The instability leads to the formation of a large number of pillars of height close to h^ax- 
As the system evolves in time, the interface self-organizes to form triangular mounds of a 
fixed slope near these pillars. These mounds then coarsen in time, with large mounds growing 
larger at the expense of small ones. In this coarsening regime, a power-law growth of W in 
time is recovered. The slope of the sides of the triangular mounds remains constant during 
this process. Finally, the system reaches a steady state with one peak and one trough and 
remains in this state for longer times. 

This behavior is illustrated in Fig^ where the interface profiles in a typical run for a 
L = 200 system starting from a flat state are shown at times t — 200 (before the onset of 
the instability), t = 4000 (after the onset of the instability, in the coarsening regime), and 
t ~ 128000 (in the final steady state). The inset shows the time-evolution of W in this run, 
as well as the results for as a function of t, averaged over 40 runs for L — 1000 samples. 
The averaged data show a power-law growth regime with /3 ~ 1/3 before the onset of the 
instability, and a second power-law coarsening regime with W tf^' , P' = 0.34 ±0.01, at long 
times. The selection of a "magic slope" during the coarsening process is clearly seen in the 
plots of FigQ] More quantitatively, the distribution of the nearest-neighbor height differences 
Si = \hi+i — hi\ is found to exhibit a pronounced peak at the selected value of the slope, and 
the position of this peak does not change during the coarsening process. The steady-state 
profile for a. L — 500 sample, also shown in Fig^ illustrates the sample-size independence 
of the results. The peak and the trough in the steady state are separated by ~ i/2 if the 
boundary condition requires the heights at the two ends of the sample to be the same. The 
occurrence of a peak and a symmetrically placed trough in the steady state is a consequence 
of using periodic boundary conditions. This symmetry is not present when other boundary 
conditions (such as "fixed" and "zero flux") are used. However, the basic phenomenology of 
mound formation, power-law coarsening and slope selection does not depend on the boundary 
conditions. 

While the saw-tooth-like surface profiles found for small c is qualitatively different from the 
self-affine morphology observed for large c, the interface width exhibits very similar power- 
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Fig. 1 - Interface profiles at three different times {t — 200, 4000, and f 28000) in a run starting from a 
fiat state for a L = 200 sample of model I with A = 4.0 and c — 0.02 (full lines). A profile for L = 500 
at t = 1.28 X 10'' for the same parameters is also shown (dashed line) with both axes scaled by 2.5. 
Double-log plots of the interface width W as a function of time t in the L = 200 run (dash-dotted 
line), and similar data averaged over 40 runs for L = 1000 samples (full line) are shown in the inset. 

Fig. 2 - Invariant height- profile for a L = 200 system with A = 4.0, c = 0.02 (full line), and a snap 
shot of the same system in steady state (dash-dotted line) . A snap shot of a L = 200 system with 
f(x) = a;/(l -I- cx), A = 4.0, c = 0.01 in the steady state (dashed line) and an invariant profile of the 
corresponding continuum equation (dotted line) are also shown. Inset: Zero-crossing of the largest 
eigenvalue «; of the stability matrix as a function of c (A = 4.0, L = 200). 



law behavior in the two cases (/? ~ /?' ~ 1/3, and a' = 1 ~ a). Thus, a measurement of 
the interface width would not distinguish between the two growth modes. A clear distinction 
between the two morphologies may be obtained from measurements of the average number 
of extrema of the height profile |14| . The steady-state profile in the mound-formation regime 
exhibits two extrema for all values of the system size L. In contrast, the number of extrema 
in the steady state in the kinetic roughening regime increases with i as a power law |14| - we 
find that for values of c for which the system is kinetically rough, e.g. for A = 4.0, c = 0.05, 
the average number of extrema in the steady state is proportional to with 8 ~ 0.83. 
This observation allows us to define an "order parameter" that is zero in the large-c, kinetic 
roughening regime and finite in the small-c, mound-formation regime. Let ai be an Ising-like 
variable, equal to the sign of the slope of the interface at site i. An extremum in the height 
profile then corresponds to a "domain wall" in the configuration of the {cTi} variables. Since 
there are two domain walls separated by '--^ i/2 in the steady state in the mound-formation 
regime, the quantity 

would be finite in the L ^ oo limit. On the other hand, m would go to zero for large L in 
the kinetically rough regime because the number of domains in the steady-state profile would 
increase with L. We find numerically that m ~ with 7 ~ 0.2 for A = 4.0, c = 0.05. 

The behavior described above may be understood from a simple stability analysis. The 
profile near the top (i = io) of a mound may be approximated as = Xq + Xi, 
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hio±j — xq — \j — l\x2, where xi is the height of the pillar at the top of the mound and X2 
is the selected slope. The condition that this profile does not change under the dynamics of 
Eq.(0J with no noise leads to the following pair of non- linear equations for xi and X2'. 

2x1 - A[l - e-"^2]/c = 0, 
3x1 - ^2 - A[l - e-'=(^i+^^)'/4]/c = 0. (6) 

These equations admit a non-trivial solution for sufficiently small c, and the resulting values 
of X\ and X2 are found to be quite close to the results obtained from numerical integration. A 
similar analysis for the profile near the bottom of a trough (this amounts to replacing X2 by 
—X2 in Eq.®) yields slightly different values for x^ and X2- The full stable profile (a fixed point 
of the dynamics without noise) with one peak and one trough may be obtained numerically 
by calculating the values of hi for which = for all z, where gi is the term multiplying 
in the right-hand side of Eq. Q . This calculation shows that the small mismatch between 
the values of X2 near the top and the bottom is accommodated by creating a few ripples near 
the top. The numerically obtained fixed-point profile for a L = 200 system with A = 4.0, 
c = 0.02 is shown in Fig|21 along with a typical steady-state profile for the same system. 
The two profiles are found to be nearly identical. Fig|21 also shows a snap shot of a L = 200 
system with j(x) — xjiX ~\- cx), A = 4.0, c = 0.01 in the steady state and an invariant profile 
of the corresponding continuum equation, obtained using the procedure of Ref. |15j . These 
results show that the steady-state properties for the two forms of jix) are very similar, and the 
continuum equation admits invariant solutions that are very similar to those of the discretized 
models. 

The local stability of the mounded fixed-point may be determined from a calculation of 
the eigenvalues of the matrix Mij — dgi/dhj evaluated at the fixed point. We find that the 
largest eigenvalue of this matrix crosses zero at c = ci(A) (see inset of Fig.©), signaling an 
instability of the mounded profile. The structure of Eq.Q implies that ci(A) oc A^. Thus, 
for < c < ci(A), the dynamics of Eq.l^J without noise admits two locally stable invariant 
profiles: a trivial, flat profile with hi the same for all z, and a non-trivial one with one mound 
and one trough. Depending on the initial state, the no-noise dynamics takes the system to one 
of these two fixed points. For example, an initial state with one pillar on a flat background 
is driven by the no-noise dynamics to the flat flxcd point if the height of the pillar is smaller 
than a critical value, and to the mounded one otherwise. We did not carry out a stability 
analysis of mounded fixed-point solutions of the continuum equation (see Fig|2Jl because doing 
such analysis without discretizing space would be extremely difficult. 

In the presence of the noise, the perfectly flat fixed point transforms to the kinetically 
rough steady state, and the non-trivial fixed point evolves to the mounded steady state shown 
in FigH A dynamical phase transition at c = C2(A) < ci(A) separates these two kinds of 
steady states. To calculate C2(A), we start a system at the mounded fixed point and follow its 
evolution according to Eq.Q) for a long time (typically t — 10^) to check whether it reaches 
a kinetically rough steady state. By repeating this procedure many times, the probability. 
Pi (A, c), of a transition to a kinetically rough state is obtained. For fixed A, Pi increases rapidly 
from to 1 as c is increased above a critical value. Typical results for Pi as a function of c for 
A = 4.0 are shown in the inset of FigO The value of c at which Pi ~ 0.5 provides an estimate 
of C2- Another estimate is obtained from a similar calculation of P2(A, c), the probability that 
a flat initial state evolves to a mounded steady state. As expected, P2 increases sharply from 
to 1 as c is decreased (see inset of FigEJ, and the value of c at which this probability is 0.5 is 
slightly lower than the value at which Pi = 0.5. This difference reflects finite-time hysteresis 
effects. The value of C2 is taken to be the average of these two estimates, and the difference 
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Fig. 3 - Critical values of the control parameter c as functions of A: ci of model I (circles), C2 of model 
I (triangles), and C2 of model II (diamonds). Inset: The probabilities Pi (circles) and P2 (inverted 
triangles) defined in text, as functions of c for model I with A = 4.0, L = 200. 

Fig. 4 - Interface profiles at three times {t = 10^, 10^, and 3.10*') in a run starting from a fiat state 
for a L = 200 sample of model II with A = 2.0 and c = 0.005 (fuU lines). A profile for L = 500 at 
t = 3.10" for the same parameters is also shown (dashed line) with both axes scaled by 2.5. Double-log 
plots of the interface width as a function of time t in the L = 200 run (dashed line), and similar 
data averaged over 40 runs for L — 1000 samples (full line) are shown in the inset. 

between the two estimates provides a measure of the uncertainty in the determination of C2. 
The phase boundary obtained this way is shown in FigEl along with the results for ci(A). 

Our numerical results for model II are very similar. Height profiles of a L = 200 sample 
with A = 2.0, c = 0.005 at three times (before the onset of the instability, during coarsening, 
and at the steady state) are shown in Fig^J along with the profile of a L = 500 sample in the 
coarsening regime. The inset shows the time-evolution of the interface width in the L — 200 
run and also the average over 40 runs for L — 1000 samples. The phase diagram for this 
model (see FigEJ is very similar to that of model I, while the coarsening exponent appears to 
have a higher value, (3' — 0.50 ± 0.01. 

The qualitative behavior found here is very similar to that in canonical first-order tran- 
sitions in two- and three-dimensional equilibrium systems. In a time-dependent Ginzburg- 
Landau description |lfij of the dynamics of a system exhibiting a first-order transition, the 
no-noise dynamics exhibits two locally stable fixed points (corresponding to the uniformly 
ordered and disordered states) in a range of temperatures near the transition temperature. In 
the presence of noise, the system selects one of the phases corresponding to these two fixed 
points, except at the transition temperature where both phases coexist. The local stability 
of the mean-field ordered and disordered states leads to finite-time hysteresis effects near the 
transition temperature. The behavior we find is very similar, with the rough and the mounded 
states corresponding to the disordered and the ordered states of equilibrium systems and c 
playing the role of the temperature. In analogy with the behavior of equilibrium systems, we 
find hysteresis and coexistence of rough and mounded morphologies near c = 02- 

In summary, we have shown that a nonlinear instability in a class of surface growth models 
leads to mound formation and power-law coarsening with slope selection via a dynamical phase 
transition. This mechanism of mound formation is very different from the conventional ES 
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mechanism. The ES instabihty is usually modeled [71|S] as a linear one arising from a V^/i 
term with a negative coefRcient. Such a term is clearly absent in our models. Consequently, a 
flat interface is locally stable in our models. Since the non-equilibrium surface current in our 
models vanishes for all values of a constant slope, the slope selection we find is a true example 
of nonlinear pattern formation. In contrast, slope selection occurs in ES-type models only 
when the surface current vanishes at a specific value of the slope jZ] . 

In view of the observation |11II17|| that spatial discretization may drastically affect the be- 
havior of continuum growth equations, it is interesting to enquire whether the truly continuum 
Eq.m with |V/i'p replaced by /(IV/i'p) would exhibit the same behavior as that found for the 
discrete models studied here. While we can not provide a rigorous answer to this question, we 
expect the continuum equation to behave in a similar way. It was shown in Ref. 11 that the 
nonlinear instability of the discretized Eq.|(2J) can not be avoided by making the integration 
time step sufficiently small, or by using more accurate left-right symmetric definitions of the 
lattice derivatives. There is no a priori reason to believe that certain left-right asymmetric 
discretization schemes 18 that avoid this instability provide a more accurate representation 
of the behavior of the continuum equation. Since the instability does not necessarily imply 
a divergence of the height ^J, the rigorous proof ^Sl of the boundedness of the solutions 
of Eq.(^ without noise does not rule out its occurrence in the continuum equation. Finally, 
the existence of mounded fixed-point solutions of the continuum equation (see Figl^J strongly 
suggests that its behavior is not qualitatively different from that of the discrete models. 

* * * 

We thank S. Das Sarma for helpful discussions and SERC, IISc for computational facilities. 
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